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Abstract 

We present results of the application of the anisotropic hydrodynamics (aHydro) framework 
to (2+l)-dimensional boost invariant systems. The necessary aHydro dynamical equations are 
derived by taking moments of the Boltzmann equation using a momentum-space anisotropic one- 
particle distribution function. We present a derivation of the necessary equations and then pro- 
ceed to numerical solutions of the resulting partial differential equations using both realistic smooth 
Glauber initial conditions and fluctuating Monte-Carlo Glauber initial conditions. For this purpose 
we have developed two numerical implementations: one which is based on straightforward inte- 
gration of the resulting partial differential equations supplemented by a two-dimensional weighted 
Lax-Friedrichs smoothing in the case of fluctuating initial conditions; and another that is based 
on the application of the Kurganov-Tadmor central scheme. For our final results we compute the 
collective flow of the matter via the lab-frame energy-momentum tensor eccentricity as a function 
of the assumed shear viscosity to entropy ratio, proper time, and impact parameter. 



PACS numbers: 12.38.Mh, 24.10.Nz, 25.75.Ld, 25.75.-q 



I. INTRODUCTION 



The goal of ultrarelativistic heavy ion collision experiments at the Relativistic Heavy Ion 
Collider at Brookhaven National Laboratory (RHIC) and the Large Hadron Collider (LHC) 
at CERN is to create a tiny volume of matter (~ 1000 fm 3 ) which has been heated to a 
temperature exceeding that necessary to create a quark-gluon plasma. Early on it was shown 
that ideal relativistic hydrodynamics is able to reproduce the soft collective flow of the matter 
and single particle spectra produced at RHIC [H-S]. Based on this there was a concerted 
effort to develop a more systematic framework for describing the soft collective motion. 
This effort resulted in a number of works dedicated to the development and application of 
relativistic viscous hydrodynamics to relativistic heavy ion collisions [oTffij]. 

One of the weakness of the traditional viscous hydrodynamics approach is that it relies 
on an implicit assumption that the system is close to thermal equilibrium which implies 
that the system is also very close to being isotropic in momentum space. However, one 
finds during the application of these methods to relativistic heavy ion collisions that this 
assumption breaks down at the earliest times after the initial impact of the two nuclei due 
to large momentum-space anisotropies in the Pt~Pl plane which can persist for many fm/c 
[27] . In addition, one finds that near the transverse and longitudinal edges of the system 
these momentum-space anisotropies are large at all times [2"Tl - f3"T] . Similar conclusions have 
been obtained in the context of strongly coupled systems where it has been shown using 
the AdS / CFT correspondence one achieves viscous hydrodynamical behavior at times when 
the system still possesses large momentum-space anisotropies and that these anisotropies 
remain large throughout the evolution [32H38] . Based on these results one is motivated to 
obtain a dynamical framework that can accommodate potentially large momentum-space 
anisotropies. 

In this paper we follow up recent work which aims to extend the applicability of space- 
time evolution models for the bulk dynamics of a quark-gluon plasma to situations in which 
there can be large momentum-space anisotropies. Initial studies along this direction focused 
on boost-invariant expansion in systems which were transversally homogeneous [39, 40J. 
The motivation and conceptual setup of Refs. (39j HQ] were similar in the sense that they 
both relaxed the assumption of the system being nearly isotropic in momentum space; how- 
ever, there was a key conceptual difference in the derivation of the resulting dynamical 
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equations. In Ref. [39J an entropy source was postulated which satisfied the minimal con- 
straints necessary in the limit of small momentum-space anisotropy and then the authors 
assumed a constant rate of isotropization regardless of the local typical momentum of the 
plasma constituents. In Ref. jlQ] the equations of motion were derived by taking moments 
of the Boltzmann equation and supplemented by a requirement that in the limit of small 
momentum-space anisotropy these equations reproduced those of 2nd-order Israel-Stewart 
viscous hydrodynamics (HH13]. The result of this matching was that the relaxation rate of 
the system was necessarily proportional to the local hard momentum scale]]] This allowed 
the authors of Ref. [40J to smoothly match onto 2nd-order viscous hydrodynamics when the 
system was nearly isotropic in momentum space. 

The phenomenological consequence of these two different results for the relaxation rate 
is quite important. If the relaxation rate is proportional to the local hard momentum 
scale, then one expects a slower relaxation to isotropy when the local hard momentum 
scale is reduced. This occurs at late times in the one-dimensional case since the local 
hard momentum scale is dynamically lowered due to expansion. Even more importantly, 
having a relaxation rate which is proportional to the hard momentum scale has important 
consequences for the evolution of the matter near the longitudinal and transverse edges of 
the system where the local temperature is also initially lower. The first demonstration of 
this effect was in Ref. [28] which studied the one dimensional non-boost invariant evolution 
of a system which was transversally homogeneous. This work followed similar developments 
in Ref. [2H] where a constant relaxation rate was assumed. A comparison of the results of 
these two papers shows that one sees much larger momentum-space anisotropics at large 
spatial rapidity being developed if one uses a relaxation rate which is proportional to the 
local hard momentum scale. 

Since these works were published, the anisotropic hydrodynamics methodology has been 
extended to include boost-invariant transverse dynamics [HI US]; however, these papers once 
again assumed a fixed rate of relaxation to isotropy. In this paper we study the effect of 
using a more realistic relaxation rate which is proportional to the hard momentum scale 
[10] . thereby allowing a smooth matching to 2nd-order viscous hydrodynamics. We present 
a derivation of the necessary equations and then proceed to numerical solutions of the 
1 In this context the hard momentum scale corresponds to the typical average momentum scale of the 

particles of the system. When one has local isotropic thermal equilibrium, the average momentum scale 

corresponds to the temperature of the system. 
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resulting partial differential equations using both realistic smooth Glauber initial conditions 
and fluctuating Monte-Carlo Glauber initial conditions. For this purpose we have developed 
two numerical implementations: one which is based on straightforward integration of the 
resulting partial differential equations supplemented by a two-dimensional weighted Lax- 
Friedrichs smoothing in the case of fluctuating initial conditions; and another that is based 
on the application of the Kurganov-Tadmor central scheme. For our final results we compute 
the collective flow of the matter via the lab-frame energy-momentum tensor eccentricity as a 
function of the assumed shear viscosity to entropy ratio, proper time, and impact parameter. 
We also present results for the dependence of the momentum-space anisotropy in the full 
transverse plane and show that in regions where the temperature is low one can develop 
sizable momentum-space anisotropies. As a control test we compare with 2nd-order viscous 
hydrodynamics in the limit of small shear viscosities and demonstrate that the aHydro 
framework is able to reproduce the temperature and flow profiles obtained from 2nd-order 
viscous hydrodynamics in this limit. 

The structure of the paper is as follows: In Sec. [XT] we introduce the tensor basis we 
will use in the case that the system is anisotropic in momentum space and derive the par- 
tial differential equations necessary for the dynamical evolution by taking moments of the 



Boltzmann equation. In Sec. Ill we present the types of smooth initial conditions we will 



use. In Sec. IV we introduce the three numerical algorithms (centered differences, weighted 
LAX, and hybrid Kurganov-Tadmor) we will we use to solve the resulting partial differen- 
tial equations. In Sec. [V] we compare with 2nd-order viscous hydrodynamics for non-central 



collisions and present our final results. In Sec. VI we present our conclusions and a future 
outlook. Finally, in three appendices we include a comparison of entropy production in 
2nd-order viscous hydrodynamics and aHydro, some numerical checks of convergence etc., 
and a brief rederivation of the 0+ld Bjorken model using our tensor formalism. 



II. KINETIC THEORY APPROACH TO ANISOTROPIC HYDRODYNAMICS 

In this section we describe our theoretical framework for describing relativistic plasmas 
which are anisotropic in momentum-space. Our setup is based on the kinetic theory approach 
to non-equilibrium systems [JT]. There are different methods for constructing approximate 
solutions of the relativistic Boltzmann equation |JT] ■ The most well-known approach is due 
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to Israel and Stewart [12J H6] . In this approach one expands the distribution function around 
a local thermal equilibrated distribution function, f eq (x,p), in terms of a series of irreducible 
Lorentz tensors q of particle momentum p^ 

f(x,p)= /eq(:C,p)(l + 0(x,t)), 

= /eq(x,p) (1 + C(x,t) + C/ J>> + W V + C^ AP <» A> + •••), (2-1) 

where the angle brackets above stand for symmetrized tensors which are orthogonal to 
the fluid four- velocity [26j HI]- The thermal equilibrium distribution function has the 
functional form 



cq 



0X1)1 T(i) J +a 



(2.2) 



where a = ±1 gives Fermi-Dirac or Bose-Einstein statistics and a = gives Maxwell- 
Boltzmann statistics. 



The distribution function (2.1) is usually expanded until second order, i.e. just keeping 
the terms 1, p^\ and p^p v \ An important aspect in the construction of irreducible tensor 
basis is the decomposition of the four-momentum p M of a particle in Minkowski space. One 
assumes the existence of a time-like normalized vector field u^(x) (which is identified with 
the fluid velocity) and an operator A^ u which is symmetric, traceless and orthogonal to 
u^(x) such that p M = Eu^ + A^ u p v [26, 41 J . This decomposition allows one to have an 
irreducible nth-rank tensor basis which is complete and orthogonal [26 : , |4"T] . 

An alternative but equivalent treatment for expanding the distribution function in terms 
of an irreducible nth-rank tensor basis was developed by Anderson (17] . This method instead 
decomposes the four-momentum p^ of a particle as 

3 

p» = Eu» + J2Pi x % > ( 2 - 3 ) 

i=i 

where is the fluid velocity and xf is a set of orthonormal vectors which are spacelike and 
orthogonal to With this decomposition one can also find a suitable irreducible tensor 
representation [47] . We will follow this decomposition closely since it is the most convenient 
2 We point out that in the original approach by Israel and Stewart, the decomposition basis is not orthogonal 
and therefore, the exact form of the transport coefficients cannot be obtained once the expansion is 
truncated. Recently, Denicol et al. showed how to correct this and expand properly the distribution 
function in terms of a complete and orthogonal set of irreducible tensors of a particle with momentum 
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vector basis for a system which is anisotropic in momentum-space along some preferred 
direction(s). 



In the rest of this section, we use the vector basis decomposition (2.3) to construct 
2nd-rank tensors. As a particular case, we construct the energy-momentum tensor for a 
(2+l)-dimensional boost invariant anisotropic plasma and derive the dynamical equations 
of motion by taking moments of the Boltzmann equation. Our discussion is restricted to 
the case of vanishing chemical potential. 



A. Vector Basis 

In this paper we will concentrate on systems which possess a preferred direction associated 
with a single direction in momentum-space. It is possible to construct a tensor basis which 
allows for multiple anisotropy directions; however, we restrict our considerations to this 
simpler case since taking into account the momentum-space anisotropy along the beamline 
direction is of particular importance for heavy-ion phenomenology. To begin, we will specify 
a tensor basis which is completely general and not subject to any symmetry constraints and 
then add the necessary symmetry constraints when needed. 

A general tensor basis can be constructed by introducing four 4-vectors which in the local 
rest frame (LRF) are 

*&rf = «£rf = (1, 0, 0, 0) 
-^t,lrf = x lrf = (0' 1' 0) 0) 

^,LRF^LRF = (0,0,1,0) 

*£lrf = *LRF = (0, 0, 0, 1) . (2.4) 

These 4-vectors are orthonormal in all frames. The vector Aq is associated with the four- 
velocity of the local rest frame and is conventionally called and one can also identify 
X± = x M , X% = y^-, and A3 = as indicated above. We will use the two different labels for 
these vectors interchangeably depending on convenience since the notation with numerical 
indices allows for more compact expressions in many cases. Note that, in the lab frame the 
three spacelike vectors X? can be written entirely in terms of Aq = w M . This is because 
Af can be obtained by a sequence of Lorentz transformations/rotations applied to the local 



rest frame expressions specified above. We will return to this issue and construct explicit 
lab-frame representations of these four- vectors later. 

Finally, we point out that one can express the metric tensor itself in terms of these 
4-vectors as 

3 

1=1 

In addition, the standard transverse projection operator which is orthogonal to Xq can be 



rewritten in terms of the vector basis (2.4) as 



A ,u = ^ _ x » x » = _ J2 X»X\ , (2.6) 

i=i 

such that u^A^ = u u A flu = 0. We note that the spacelike components of the tensor basis 
are eigenfunctions of this operator, i.e. X^A^ = X\ . 

B. 2nd-rank Tensors 

A general rank two tensor can be decomposed using the 4-vectors Xg. In general there 
are sixteen possible terms 

3 

^(t, X )= c ^ x " x h 

a,/3=0 

3 3 
= CooX^Xq + CjjXfXf + c apX£Xp , 

1=1 a,0=O 

3 3 

= W" + (<* + c oo) X i X i + E °^ X a X p > (2-7) 

1=1 «,/3=0 
= "ii 

where it is understood that the coefficients c a p now contain all of the space-time dependence. 



C. 2nd-rank symmetric Tensors 

If a two tensor is symmetric under the interchange of /i and v then c a/ 3 = cp a and we can 
write 

3 3 

A^(t, x) = c 00 <r + £ W + £ c a/3 (X£X£ + X£X^) . (2.8) 

i=l q,/3=0 
o>/9 

and there are only then ten independent terms. 
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1. Energy-Momentum Tensor for Ideal Hydrodynamics 



Since the energy-momentum tensor is a symmetric tensor of 2nd-rank, Eq. (2.8) can be 
used 

3 3 

T^(t, x) = fa^ + £ ^ + £ + X*X$ , (2.9) 

1=1 <*,/3=0 

where we have relabeled the coefficients for this purpose. In the local rest frame we can 



identify the basis vectors via (2.4) and we have that = £ and T"lrf = ^ where £ is the 
energy density and V% is the pressure in z-direction and all other components vanish. If the 



system is locally isotropic as is the case for ideal hydrodynamics then V% = V. From (2.9) 
we have T t °Stt = £ = t nn and T^ RF = V = —too + U% an d since all off-diagonal components 



vanish we have t Qj g = for all a ^ (3. This allows us to write 

3 

T^(t, x) = EgT + (V + £)J2 X K , 

i=l 

= £g>"+(V + £)(X$Xg- ! r) 

= (e + v)x*x%-v<r, 



(2.10) 



where in going from the first to second line we have used Eq. (2.5). Using the conventional 



notation that Xq = u M we obtain 



= (£ + V)u"u u - Vg 



(2.11) 



in agreement with the expected result. For later use we also note that 



T% = T = £ - 3P . 



(2.12) 



2. Energy-Momentum Tensor for A zimuthally- Symmetric Anisotropic Hydrodynamics 

In the bulk of this paper we will consider systems for which the momentum-space particle 
distribution is azimuthally symmetric while the rotational symmetry in the p±~Pl plane is 
broken. From here on we will refer to this as "azimuthally-symmetric" which only implies 
an assumed symmetry in momentum-space and not in configuration space. In the case of 
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azimuthally-symmetric anisotropic hydrodynamics we have 



00 c _ ± 

LRF — c — H)0 , 



J LRF — ' - 



TVV — T> 
LRF — ' 



LRF 



—too + ^11 > 
— ioO + ^22 j 
*Pl = —too + ^33 , 



(2.13) 



and due to the azimuthal symmetry in momentum-space we must have in = i22 which gives 
four equations for our four unknowns. Solving for the coefficients i one obtains 

2 

r^(i, x) = £$r + (v± +£)J2 x i x i + ( Vl + £ ) x s x z » 

i=i 

3 



£<T + (P± + £) X i X i + - P ±) X S X 3 



(£ + v±)x»xx - VlsT + {v L - P±)x%xy 



(2.14) 



Relabeling Xq = m m and X% = to agree more closely with the notation of Ref. [IS] we 
obtain 

Tf »> = + r ± ) u ^ _ <p ±g ^ + (p L _ Va_)z»z v , (2.15) 



which in the limit that V± = Vl = V reduces to (2.11). We again note for later use that 



T% = T = £ - 2V ± - V L 



(2.16) 



D. Explicit Forms of the Basis Vectors 

In the lab frame the three spacelike vectors X? can be written entirely in terms of Xq = 
u u '. This is because Xf can be obtained by a sequence of Lorentz transformations/rotations 
applied to the local rest frame expressions specified above. To go from the lab frame to LRF 
we can apply a boost along the z-axis followed by a rotation around the z-axis and finally a 
boost along the x-axis, i.e. mlrf — L x (ip)R z (6)L z (i9)u [IB] . This specific transformation is 
chosen in order to ensure that the four-vector z^ has no transverse components in all frames. 



To find the necessary vectors in the lab frame based on the LRF expressions (2.4) we apply 
the inverse operation X£ )LAB = (L X R Z L Z )~ X X^ LBF = {L z )~ x {R z )~ x {L x )~ x X^j JRF which is 



explicitly given by 



a, LAB ■ 



/ coshtf sinh$ \ / 1 

1 cos0 

1 sin0 

\ smli?? cosh$ / \ 

v " 



\ / cosh0> sinh0> \ 

- sin sinh ijj cosh ip 

cos0 1 

1/ \ 1 / 

V 



which gives 



m° = cosh -0 cosh $ . 
u 1 = sinh ^! cos , 
m 2 = sinh -0 sin , 
u 3 = cosh ip sinh $ , 

y° = o, 

y 1 = - sin , 
?/ 2 = cos , 

y 3 = o, 



x° = sinh -0 cosh $ , 
x 1 = cosh ip cos , 
x 2 = cosh -0 sin , 
s 3 = sinh ip sinh , 

z° = sinh d , 

= 0, 
z 2 = 0, 
z 3 = cosh i? . 



a,LRF • 



(2.17) 



(2.18) 



In the limit that the system is boost invariant one can identify $ = <r, where £ is the 
spatial rapidity defined through 



t — t cosh £ . 
2 = r sinh q , 



(2.19) 



where r = \/t 2 — z 2 is the proper time. In the remainder of the paper when we refer to a 
boost-invariant system we will use r and <r as the longitudinal coordinates. 



E. Dynamical Equations 

In this section we derive the dynamical equations of motion by taking moments of the 
Boltzmann equation jH] 

I?d ll f(x,p) = -C\f\. (2.20) 
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The moments are defined by multiplying the left and right hand sides of the Boltzmann 
equation by various powers of the four-momentum and then averaging in momentum space. 
This can be achieved via the n th moment integral operator 



l n = j dxp^P 112 ■■■P fln , (2.21) 

where n > is an integer and 

fo s Swf ,( W " m2) 26{f) = lw?¥»- (222) 

F. Zeroth moment of the Boltzmann equation 

The zeroth moment of the Boltzmann equation results from applying X to both sides of 



(2.20) 



dx P^d^f = J , 

f d 3 p p» 

d »Jj2^p^ j - Jo > 

<V = J , (2.23) 

where J n = —X n C[f}. Note that we can rewrite the left hand side of the last expression as 
j M = nu^ where n is the particle number density in the local rest frame. Expanding we find 

d^f = Dn + n9, (2.24) 

where 

D EE , 

0EE3X' (2.25) 

allowing us to write a general expression for the zeroth moment of the Boltzmann equation 

Dn + n6 = J Q . (2.26) 

G. First moment of the Boltzmann equation 

The first moment of the Boltzmann equation is equivalent to the requirement of energy 
and momentum conservation [H] 

d^T^ = , (2.27) 
11 



where T^ v is the energy momentum tensor. In the following we derive evolution equations 
under different assumptions about the degree of symmetry of T^ v . 



1. Ideal hydrodynamics 



To begin we use the general form of the energy-momentum tensor for an isotropic system 



given in Eq. (2.11) to obtain 



dpT 1 *" = u v D{8 + V)+ u v {8 + V)6 + {8 + V)Du v - d v V 



(2.28) 



where D and 6 are defined in Eq. (2.25). 



Canonically one takes projections of d^T^ = parallel and perpendicular to u M . The 
parallel projection is obtained via u^d^T^ which gives 



UudpT^ = D(8 + V) + (8 + V)6 + (8 + V)u v Du v - DV = 
= D8 + (8 + V)0 = 



(2.29) 



where we have used u v u v = 1 and u u Du u = ^D(u v u v ) = 0. This gives us our first equation 



for ideal hydrodynamics. For the transverse projection we use A MW defined in Eq. (2.6) which 
satisfies A au u u = 0. This gives 



A a u d fl T fJ ' v = (8 + V)A a v Du v - A a u d u V = . 



(2.30) 



Using the explicit form for A a v = g a u — u a u v one obtains A a u Du u = Du a . We can addi- 
tionally define 

3 

V Q = A a u d v = -J2 X«X u pd u , (2.31) 



which is the gradient in the spacelike directions. Putting this together with Eq. (2.29) one 
obtains the following two equations 



D8 + (8 + V)6 = , 
(8 + V)Du a - V a V = . 



(2.32) 



In the second case a should be a spacelike index such that we have four equations in total 
which should be supplemented by the equation of state which can be expressed in the form 
of a constraint on the trace of the energy momentum tensor T M M = T = 8 — 3V. 
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2. Ideal Boost Invariant Dynamics with Transverse Expansion 



In this section we briefly review what happens when the system is boost invariant and 
we allow for inhomogeneities and flow in the transverse direction. In this case we have from 



(2.18) 



ur 



(cosh if) cosh q, sinh vb cos 0, sinh if> sin 0, cosh if) sinh q) . 



(2.33) 



It is convenient at this point to relabel the components of "uF" as 



m m = (uq cosIk, u x , u y , uq sinh^) 



(2.34) 



where the constraint Uq = 1 + u x + should be satisfied. Changing to proper time and 
spatial rapidity we obtain u T = u , u q = 0, and we have 



D = u^d^ = u d T + u ± ■ V_l , 
9 = d.v? = d T u + V±- Ul + 



For the transverse gradient it is convenient to rewrite 



Uq 
T 



(2.35) 



V i = A l u d" = {g\ - v}u u )d v = d i - u l D 



(2.36) 



such that the second equation in (2.32) can be expanded into three equations 



which together with 



{£ + V)Du x + u x DV + d x V = 
(£ + V)Du y + u y DV + d y V = 
(8 + V)Du + u DV -d T V = 

D£ + (£ + v)e = a, 



(2.37) 



(2.38) 



would seem to give four equations for our four unknowns (£, V, u x , and u y since ul 



1 + u x + ul); however, upon inspection one finds that Eqs. (2.37) are not independent 



since uo times the third equation is equal to u x times the first plus u y times the second. 
We, therefore, have a choice of which equations to use and one can pick two of the three 



equations from (2.37), e.g. the first two. The final equation is then provided canonically by 



the equation of state which specifies, e.g., the energy density as a function of the pressure. 
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H. Azimuthally-Symmetric Anisotropic Hydrodynamics 



We now proceed to the derivation of the dynamical equation for azimuthally-symmetric 
anisotropic hydrodynamics. We remind the reader "azimuthally-symmetric" means that 
the momentum-space particle distribution is azimuthally symmetric while the rotational 
symmetry in the p±-Pl plane is broken. To begin we use the general form of the energy- 



momentum tensor for an azimuthally-symmetric anisotropic system given in Eq. (2.15) to 
obtain 

d^ v = u v D{£ + V±) + u u (£ + V ± )9 + (8 + Vi_)Du v - d u Vx 

+z»D L (V L - V ± ) + z"(V L - V ± )6 L + (V L - V ± )D L z» = , (2.39) 

where 

9 L = d^. (2.40) 

As before we take projections of d^T^ v = parallel and perpendicular to u M . The parallel 
projection is obtained via u u d /J T flv which gives 

u v d^ v = DS + (£ + V±)6 + (V L - V ± )u u D L z» = , (2.41) 

where we have used u v u v = 1, u v Du v = \D{u u u u ) = 0, and u v z u = 0. This gives us our 
first equation for azimuthally-symmetric anisotropic hydrodynamics. 



For the transverse projection we use A^ u defined in Eq. (2.6) which satisfies A au u u = 
and A av z v = z a . This gives 

A a v d^ v = (S + V ± )Du a - V a V ± + z a D L (V L - V±) + z a (V L - V ± )9 L 

+ (V L - V ± )D L z a - (V L ~ V ± )u a u,D L z u = . (2.42) 

1. Boost Invariant Dynamics with Transverse Expansion 

In this case we have z T = and z v = 1/r such that 

D L = z»d„ = ^ , 

T 

9 L = d^ = 0. (2.43) 
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From the first line above we find u v Dlz u = uq/t. This allows us to simplify the parallel 
projection to 

D£ + (£ + V ± )6 + (V L -V ± )— = 0. (2.44) 

r 



The transverse projections can also be simplified to 



{£ + V ± )Du a + u a DV ± + d a V ± + (V L - Vx) I - -^u a J = , (2.45) 
from which we can then obtain three equations 

(8 + V ± )Du x + u x DV± + d x V ± + (V ± -V L )^ = 0, 

T 

(S + V ± )Du y + u y DV ± + d y Vx + (V ± -V L )^ = 0, 
(£ + V ± )Du + u DV ± - d T V ± + (V± - V L ) -± = . (2.46) 

T 

As was the case with ideal hydrodynamics, we see that Uq times the third equation is equal 
to u x times the first plus u y times the second so that it is redundant. This leaves us with 
the following three equations 

D£ + {£ + V±)9 + (V L - V 
(£ + V±)Du x + 8 X V± + u x DV± + (V± - V L 

(£ + V ± )Du y + d y V ± + u y DV ± + (V ± - V L ) ^ = . (2.1 7) 

I. Distribution function for azimuthally-symmetric systems 

We next consider the one-particle distribution function / in the local rest frame and show 
that in the case of a system that is locally azimuthally-symmetric in momentum space that 
it suffices to introduce one anisotropy parameter £ and a single scale A [25]. To begin we 
consider the general form 





= o, 






U U X 


= 0, 




T 


U Uy 


= 0. 


T 





/(t,x,p) = fU^P^H^Pu) ■ (2.48) 

S M1/ (t, x) is a symmetric tensor, /j SO is an arbitrary isotropic distribution function, and p M = 
p^/A, where A(t,x) is a momentum scale that can depend on space and time (the so-called 
hard momentum scale). In the case where the system is in thermal equilibrium, then /; so 
would be given by a Bose-Einstein or Fermi-Dirac distribution function. Note that the 

15 



argument of the square root in /i SO should remain greater than or equal to zero in order for 
/ to be a single- valued real function. 

If is a symmetric tensor and is diagonal in the local rest frame, we have 



CaouW + ^cuXfXZ, (2.49) 



1=1 

and if, additionally, the system is symmetric under x •f-)- y then en = C22 = c±± and we have 

2 



1=1 



c o«'V - c ±± A^ + (c 33 - c±x)XgXZ . (2.50) 



Using our ability to redefine A — > y / cooA in Eq. (2.48) we can rescale our coefficients. 
Defining c_lj_/coo = $ and (C33 — c_lj_)/coo = a we can write compactly 

S"" = - + q// . (2.51) 

Contracting with four-momenta on both sides we find 

V^ v Vu =P 2 o + $P 2 + apl , 

= m 2 + (1 + $)p 2 + ap 2 z , (2.52) 

where we have used p$ = p 2 + m 2 . If we have a system of massless particles then 

Plt S^p v = (1 + $K + (1 + $ + a)p 2 , (2.53) 

and in this case we can once again use our ability to rescale A — > a/(1 + $)A and defining 
1 + f = (1 + $ + a)/(l + $) we obtain 

p M H^=pi + (l + 0pL (2.54) 

which has the form of the argument of the original one-dimensional Romatschke-Strickland 
(RS) distribution function [49J. 



J. Number density and Energy-Momentum Tensor with the RS distribution func- 
tion 

Based on the results of the last section, the functional form of the RS distribution function 
for a locally azimuthally-symmetric expanding anisotropic plasma is 



/(x, P ,r) = /bs(p,& A) = / iso (V[pi + (1 + £M]/A 2 ) , (2.55) 
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where it is understood that on the right hand side £ and A can depend on space and time. 
Using this distribution function the number density is given by 



n&A)= I W y 3 f™ = 7m - (2-56) 



where rii SO (A) is the number density one obtains in the isotropic limit. 
One can also evaluate the energy-momentum tensor in the LRF 



By using the RS form (2.55) one gets the explicit components of the energy-momentum 
tensor [51] 

— T TT — £;so(A) , (2.58a) 

V ± (A,0 = \ (T xx + V y ) = 7^(£)7\ (A) , (2.58b) 

V L (A, = -T? = ^ L (0^iso(A) , (2.58c) 

where Pi SO (A) and £i SO (A) are the isotropic pressure and energy density, respectively, and 

n(0 = Ut^ + aTCt Z^ ) , (2.59a) 



2 vi + e 

R±(£) s i , (2,9b) 

Ri(£) s ? f E + Wa-n ■ (2.590 



{ V f + i 



The equation of state can be imposed as a relationship between £ iso and Pi S0 . In what 
follows we will assume an ideal equation of state which is appropriate for a conformal massless 
gas, i.e. £ iso = 3"P iso - 

K. Relaxation time approximation 

As mentioned in previous sections the dynamical equations necessary can be obtained 
by taking moments of the Boltzmann equation p^d^f = — C[f]. Here we use the relaxation 
time approximation with relaxation rate V 

C[f RS ]=p^T [/ RS (p,£,A,^)-/ cq (|p|,T)] , (2.60) 
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where <j is the spatial rapidity and we fix T such that the 2nd-order viscous hydrodynamical 
equations are reproduced in the one-dimensional transversally symmetric case [3D]. This 
requires that 



7V 



which for an ideal equation of state results in 

_ 2T(r) _ 2ft 1 / 4 (QA 
577 5?7 



(2.62) 



where fj = rj/S with 77 being the shear viscosity and 5 being the entropy density. We note 
that one could perform a matching to 2nd-order viscous hydrodynamics including transverse 
dynamics, but we have not attempted to do so. Instead we use the Id matching above and 
in the results section we show that numerical results from viscous hydrodynamics codes 
which include transverse dynamics are reproduced for small fj. That being said, we have 
no reason to expect that the linearized equations would not reproduce 2nd-order viscous 
hydrodynamics; however, this remains to be proven. 



L. Dynamical Equations of Motion 

Based on the results of the previous sections, we can derive the explicit form of the 
dynamical equations of motion for a (2+l)-dimensional boost invariant system. 



1. Zeroth moment of the Boltzmann Equation 



For the RS form the Oth moment of the Boltzmann equation (2.26) is written as 



Y^ D Z ~ 6-D(log A) - 29 = 2T (l - H 3 ^) VT+T) 



(2.63) 



where we used explicitly the functional form of particle density n (2.56) and the scattering 



kernel for relaxation time approximation (2.60) 



18 



2. First moment of the Boltzmann Equation 



Using the RS form one finds the following three equations by requiring energy-momentum 
conservation 

K'{Z)DZ + 4ft(0D(log A) = - (k($ + ^±(0) a± - (^(0 + j^MC)) 7 > 

+ ^±(0] ^ = "«± \K'A0D£ + 4K ± (Z)D(\og A) + ^(Kl(0 - 7e L (0)l , 

L r J 

M J [3^(0 + K ± (0\ D (^j = K' ± (Z)D ± Z + 47^(£)£>±(logA) , (2.64) 



where 



A_l = d T U + V_l • U ± 



) 

2 



D = u d T + -%-u± ■ V± , 
u\ 

D± = z ■ (u ± x V T ) = Mx9 9 - u y d x , (2.65) 
uj_ = (u x , u y ), and Mq = 1 + u\. 

III. INITIAL CONDITIONS 

We consider collisions of symmetric nuclei, each containing A nucleons. We will study 
both participant and binary collision type initial conditions jS2] using a Woods-Saxon distri- 
bution for each nuclei's transverse profile [53]. For an individual nucleus we take the density 
to be 

n ^ r ) = 1 + J r °-fl)/ d » 

where uq = 0.17 fm~ 3 is the central nucleon density, R = (1.12V4 1 / 3 — 0.86A -1 / 3 ) fm is 
the nuclear radius, and d = 0.54 fm is the "skin depth". The density is normalized such 
that lim^oo J d 3 rnA(i r ) = A, where A is the total number of nucleons in the nucleus. The 
normalization condition fixes n to the value specified above. From the nucleon density we 
first construct the thickness function in the standard way by integrating over the longitudinal 
direction, i.e. 

/oo 
dzn A (^x 2 + y 2 + z 2 ). (3.2) 
-oo 

With this in hand we can construct the overlap density between two nuclei whose centers are 
separated by an impact parameter vector b which we choose to point along the x direction, 

19 



i.e. 6 = bx. We choose to locate the origin of our coordinate system to lie halfway between 
the center of the two nuclei such that the overlap density can be written as 



n AB (x, y, b) = T A (x + 6/2, y)T B (x - 6/2, y) . 



Another quantity of interest is the participant density which is given by 



(3.3) 



n pa ,rt(x,y,b) = T A (x + 6/2, y) 



(r NN T B (x - 6/2, y) 



B 



B 



+ T B (x-b/2,y) 



(T NN T A (x + 6/2,|/) 
A 



• (3-4) 



For LHC collisions at ^snn = 2.76 TeV we use u^n = 62 mb and for RHIC collisions at 
a/ snn = 200 GeV we use ctnn = 42 mb. From the participant density we construct our first 
possible initial condition for the transverse energy density profile at central rapidity 

n paIt (x,y, b) 



^part g 

° " ° n part (0,0,0) ' 



(3.5) 



where Sq is the central energy density obtained in a central collision between the two nuclei. 

As an alternative initial condition for energy density one could use the number of binary 
collisions which is defined as 



n co ii(z, y, b) = a NN n AB (x, y, 6) . 

from which we obtain the binary collision energy scaling 

^ con = £ n co n(x } y } b) = n AB (x,y,b) 
n coI1 (0, 0, 0) °n AB (0,0,0) 



(3.6) 



(3.7) 



IV. NUMERICAL METHODS 



We consider both smooth and fluctuating initial conditions using three numerical algo- 
rithms. In the following two subsections we describe the implementation of each algorithm. 
In each case detailed below the code is implemented using the C programming language. 



A. Centered Differences Algorithm 



In the first algorithm which we will refer to as the "centered-differences algorithm" we 
solve Eqs. (2.63) and Eqs. (2.64) by first analytically solving for the individual proper-time 
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derivatives of the four dynamical variables: £, A, u x , and u y using Mathematica [51]. We 
then had Mathematica output, in C format, the necessary right hand sides of the four 
update equations. We then discretize space on a regular square lattice with lattice spacing, 
Ax = a. For the spatial derivatives we use centered differences except on the edges of the 
lattice where we apply either a left- or right-handed first order derivative. For the temporal 
updates we use fourth-order Runge-Kutta (RK4) with a step size of At = e. 

For smooth initial conditions the previous method suffices; however, for fluctuating initial 
conditions one finds that using centered differences introduces spurious oscillations in regions 
where there are large gradients. In order to damp these oscillations one could attempt to 
use a two-dimensional Lax-Friedrichs (LAX) update [551 EE]. In practice this amounts to 
replacing the current value of a given dynamical variable by a local spatial average over 
neighboring sites and using this as a stand in for the current value of the variable, e.g. 



£lax(t, x, y) = [f (r, x + a, y) + f (r, x-a,y) + £(r, x, y + a) + f (r, y - a)] /4, (4.1) 
and now the £ update for a temporal step of size e becomes schematically 



where RHS^ stands for the (rather complicated) right hand size of the £ update equation. 
However, such a scheme results in too much numerical dissipation. An alternative is to 
realize that the source of the spurious oscillations is the weak coupling between odd- and 
even-number lattice sites. The full LAX scheme above maximally couples these interleaving 
lattices; however, this need not be done. Instead one can weight the LAX-smoothed values 
with a weight A and combine this with the current value of the variable in question, e.g. 



The smaller the value of A, the less the numerical viscosity. In practice, we have found 
that for the aHydro equations one should take A > 0.02 in order to achieve numerical 
stability. In the results section below we use A = 0.05 which represents a factor of twenty 
decrease in the dissipation induced by LAX-smoothing. Note that, when activated, wLAX 
smoothing is implemented for all dynamical variables (A, £, u x , and u y ) after each full time 
step of e and not within each RK4 substep. We will only need to use the wLAX method for 
fluctuating initial conditions; however, in App. |B]we present numerical tests using it in the 



£(t + e, x, y) = £lax(t, x, y) + e RHS 5 (r, x, y) , 



(4.2) 



£wLAx(r, x, y) = \£lax(t, x, y) + (1 - A)f (r, x, y) . 



(4.3) 
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smooth initial conditions case in order to show that the amount of numerical viscosity in 
the wLAX case is not numerically significant. That being said, one would also like to have 
another method for handling the spurious oscillations caused by using higher-order centered 
differences. This has motivated us to also implement the Kurganov-Tadmor central scheme 
which we describe in the next subsection. 

B. The MUSCL Algorithm 

As mentioned above, when there are large gradients present in a hyperbolic partial differ- 
ential equation, the application of straightforward centered-differences scheme can lead to 
spurious oscillations. For smooth initial conditions and finite shear viscosity this is not an 
issue; however, for fluctuating initial conditions one needs a way to handle shocks and dis- 
continuities. One way to proceed is to implement the LAX method as described previously; 
however, the LAX method introduces numerical viscosity into the algorithm which scales like 
the (Ax) 2 / At so that it is not possible to take the temporal step size to zero without having 
extremely small lattice spacing to reduce the numerical viscosity. As discussed above one 
can reduce the amount of numerical viscosity by instead using the weighted LAX (wLAX) 
prescription described above; however, it is desirable to have an alternative algorithm in 
order to be sure of the results. 

For this purpose we have also implemented a "Monotone Upstream-Centered Schemes for 
Conservation Laws" (MUSCL) scheme derived by Kurganov and Tadmor [57] which has been 
extended to include nonlinear sources [HE]- This method is particularly appealing because 
it can be shown that, although it does induce some numerical viscosity, the magnitude of 
the numerical viscosity induced scales like as a power of the lattice spacing with no power of 
the temporal step size in the denominator allowing one to take extremely small time steps 
without inducing large artificial numerical viscosity. Our implementation closely follows that 
introduced by Schenke et al. [59] to solve three-dimensional relativistic ideal hydrodynamics 
equations. They have also extended the method to 2nd-order three-dimensional relativistic 
viscous hydrodynamics [201 EI] with fluctuating initial conditions. 

To explain the algorithm let us consider the simpler case of a one dimensional system of 
hyperbolic partial differential equations which can be cast into "conservative" form, i.e. 




(4.4) 
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where u is, in general, an n-dimensional vector, F is a so-called flux variable or flux function, 
and F x (u) = d x F(u). For example, if one were solving the advection equation, d t u + d x u = 
then we would have F = u and if one were solving Burgers' equation d t u + ud x u = this 
can be written in conservative form as dtu + d x (u 2 /2) = so that, in this case, F = u 2 /2. 



Given a partial differential equation of the form (4.4) Kurganov and Tadmor derived the 
following semi-discrete update equation 



du-i 



Hj+i/ 2 (t) - Hj_ 1/2 (t) 



dt Ax 
where the numerical flux function H is given by 



(4.5) 



Hj+i/ 2 {t) 



Fluj +1/2 (t) 



l j+l/2 



(t) 



3+ 



2 2 

with o| +1 / 2 (t) being the local propagation velocity in the x-direction which is given by the 
maximum of the left and right half-site extrapolated spectral radius of dF/du which is 
defined as p 



4+i/ 2 (*) = max {p(j[ (4+1/2^))) > p(jj^ 



u 



i+1/2 



(t) 



\ du 

and finally, the half-site extrapolated intermediate values u^ +1 j 2 are given by 



(4.7) 



U j+l/2 



U j+ i(t) 



Ax 



[u 



x)j+l 



(t) 



Ax 



u j+ i/2 = u j(t) + -irMj(t) ■ 



For the derivatives, u x , appearing in (4.8) one should use a total variation diminishing 



"flux-limiter" so that spurious oscillators are avoided [SD]- We follow the original paper of 
Kurganov and Tadmor and use the three-argument minmod flux-limiter |61j : 



(u 



minmodf Q ^tl , «£±^t* U ^ 
V Ax ' 2Ax 



u, 



Ax 



, 1 < 9 < 2, 



(4.9) 



where 



minmod(xi, x 2 , 



mmj{xj}, if Xj > V j 

maxjixj}, if Xj < Vj (4.10) 
otherwise 

The value of 9 controls the dissipation of the flux limiter with 6 = 1 being the most dissipative 
and 9 = 2 being the least. In this paper we follow [59] and use 9 = 1.1. For details of 
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the derivation of the Kurganov-Tadmor scheme we refer the reader to their original paper 
[oTj . As mentioned above one can extend the Kurganov-Tadmor scheme to accommodate 
nonlinear time-dependent sources. Including the possibility of a time-dependent source 
changes our one-dimensional example to 

d t u + F x (u) = J(t,u), (4.11) 



where J is a source term. Naidoo and Baboolal [SB] demonstrated that, in this case, only a 
simple modification of adding the source on the right hand side was necessary 

^ = _ H i+ i/2(t) - Hi_ 1/2 (t) + J{t ^ (4l2) 

We note that to extend the method described thus far to multiple dimensions one introduces 
flux functions for each direction, e.g. F y and F z , and includes these in the update rule by 



defining new numerical flux functions (4.6) and propagation velocities (4.7) accordingly. 



1. Applying MUSCL to aHydro 

In the case of aHydro all of the evolution equations stem from conservative systems 
with sources, therefore we can apply the general method just described. For this purpose 
we need the first and second moments of the Boltzmann equation with the RS form for the 
one-particle distribution function. The zeroth moment can be written in a conservative form 
with sources in r-q coordinates as follows 

d T f + V ± -3 ± = - J - + Jo, (4.13) 

T 



where j M = n is the particle four-current and 

J = rn iso (A) ' 



(4.14) 



is the zeroth-moment of the right-hand side of the Boltzmann equation in the relaxation 
time approximation used herein. The remaining three update equations necessary can be 
obtained from energy-momentum conservation, d^T^ = 0, giving 

d T T TT + d x T TX + d y T Ty = — \T TT + r 2 T«l , (4.15) 

T 

i-prx 

d T T TX + d x T xx + d y T xy = , (4.16) 

r 

d T T ry + d x T xy + d y T yy = _'£1_ ( 417 ) 
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Once the dynamical variables j T , T TT , T TX , T ry are updated via these equations, they can 
then can be used to construct the remaining components of j M and T^ v . In our case it is 
necessary to solve two simultaneous nonlinear equations for £ and A which will then allow 
us to determine the rest of the information necessary to proceed with the solution. To see 



how this works in practice, we first use (2.15) to write the non-vanishing components of 
and j M = nu^ explicitly 



rpTT 


= (£ + V ± 


)u°u° - V± , 


(4.18a) 


r-pri 


= (£ + V ± 


)u°u l , 


(4.18b) 


rpij 


= (S + V ± 


)u l u j , 


(4.18c) 


rpii 


= (S + V ± 


) M V + V± , 


(4.18d) 




= V L /r 2 , 




(4.18e) 


f 


= nu° , 




(4.18f) 


f 


= nu 1 , 




(4.18g) 



where % G {x, y}. Using these equations and the normalization condition = 1 + u 2 x + v? y 
one finds two nonlinear equations, similar to those obtained in Ref. 

nrT (t tx ) 2 + (T T yf 



and 



J 



n(A,0 



T TT + V±(A,0 



(4.19) 



(4.20) 



From these two equations one can numerically solve for A and £. These values can then be 
used to determine u T and u l via 



u 



u 



J' 



n(A,()T' 



(4.21a) 
(4.21b) 



Once determined, these components of the four- velocity together with the values of A and 



£ can be used to determine all remaining variables in (4.18). 

The only remaining ingredient necessary for the Kurganov-Tadmor algorithm to be im- 
plemented fully is to determine the local propagation velocities «j +1/ / 2 (^)- These are obtained 
by evaluating the eigenvalues of the 4x4 Jacobian of j T , T TT , T TX , T Ty . As was the case in 
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Ref. |59j, with some work and a little bit of help from Mathematica, one finds that two 
of the four eigenvalues are degenerate and equal to u l ju T and the other two are given by 



A, 



± 



A± VB 
D 



(4.22) 



with 



A 
B 
D 



uV(l - 

,2 



,2 



' 2 2 

ut — uf 



(«? - l)v 2 , 



and 



dVi 



+ 



n 



88 £ + V ± dn ' 
Using an ideal equation of state for which £ iso = 3"Pi SO one obtains 

2 = 1 2^(0 + 3(1 + 0^1(0 4(i + o ^(0^(0-^(0^1(0 

V UJ 3 2^(0 + 3(1 + 0^(0 3^(0+^ ± (0 2^(0 + 3(1 + 0^(0 
In this function both terms individually diverge in the limit that £—>•(), however, these diver- 
gences cancel to give a finite result of lim^of 2 = 2/5. It has other limits of limg_;._i v 2 



(4.23a) 
(4.23b) 
(4.23c) 

(4.24) 
(4.25) 



and lini£_ 





1/2. Using the now known eigenvalues one finds that the maximum value 



of the four eigenvalues is given by 



p = |max(Aj)| 



\A\ 



B 



D 



(4.26) 



Using the above scheme one can evolve the aHydro system with fluctuating initial 
conditions; however, there is a caveat, namely that the linearly interpolated intermediate 
values of j T , T TT , T TX , and T Ty determined via (4.8) may not have real-valued solutions 



for A and £ using Eqs. (4.19) and (4.20). In practice, we find that it is necessary to use 



extremely fine lattices in order to ameliorate this problem. Alternatively, we have found that 
instead of extrapolating the four variables j T , T TT , T TX , and T Ty to the half-sites, one can 
instead extrapolate the current values of A and £ to the half-sites for use in evaluating the 
flux functions. In addition, we have found that in practice it is necessary to use a "hybrid" 
algorithm in which the centered-differences scheme described in the previous subsection is 



used as the initial guess for the nonlinear root finder which solves Eqs. (4.19) and (4.20). 



This is necessary, in particular, in regions where £ ~ since Eqs. (4.19) and (4.20) have 
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two solutions which become very close together causing the root finder to oscillate between 
the two solutions. The predicted value from the centered-differences scheme predicts for the 
nonlinear root finder which solution to use in this case. We will refer to this method as 
"Hybrid Kurganov-Tadmor" . 



V. RESULTS 



In this section we present results for the time evolution of the matter generated in heavy 



ion collisions at LHC energies using the aHydro evolution equations (2.63) and (2.64). 
For the results presented here we assume a ideal gas of quarks and gluons with Nf = 2 so 
that there are A^jof = 37 degrees of freedom. For our numerical tests and results we will 
concentrate on the spatial and momentum-space ellipticities, e x 

* - <y 2 -* 2 >£ , (5.i) 



<x 2 + y 2 >£ 



and e p is defined in the lab frame via 



e„ = — , (5.2) 

where <x 2 >£ and <y 2 >£ are the proper-time dependent average values of x 2 and y 2 weighted 
by the energy density 

<x 2 > £ = J\f [ x 2 £(r,x,y), (5.3) 

Jx,y 

and the averages in the momentum-space ellipticity represent unweighted integrals over the 
transverse directions. 

Note that the normalization M is arbitrary since it cancels in the ratio we are computing. 
These definitions are the conventional ones from the literature j62] which, unfortunately, are 
slightly inconsistent since e x is defined in the local rest frame and e p in the lab frame. It 
would be more consistent to weight the spatial average by T TT ; however, to be consistent 
with the existing literature we will use the definition weighted with the energy density in 
the local rest frame. 

We concentrate on the ellipticities since, as we will see, large momentum-space anisotropies 
are developed during the evolution of the system. Such large momentum-space anisotropies 
cast doubt on the naive application of Cooper-Frye [63] and linearly-corrected Cooper- 
Frye [S3]. We, therefore, postpone the implementation of freeze out until we can allow 
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FIG. 1. (Color online) Comparison of aHydro isotropic temperature and flow profiles with 2nd- 
order viscous hydrodynamics code for Airrj/S = 0.1 and 6 = 7 fm. Lattice size used was 109 x 109 
with a = 0.394 fm, e = 0.01 fm/c, t = 0.25 fm/c, A = 600 MeV, and Co = 0. For the transverse 
profile Glauber binary collision scaling was used. 



4;tr|/S = 10 47ir|/S = 10 




y [fm] y [fm] 



FIG. 2. (Color online) Comparison of aHydro isotropic temperature and flow profiles with 2nd- 
order viscous hydrodynamics code for 4irri/S = 10 and 6 = 7 fm. Lattice size used was 109 x 109 
with a = 0.394 fm, e = 0.01 fm/c, t = 0.25 fm/c, A = 600 MeV, and £ = 0. For the transverse 
profile Glauber binary collision scaling was used. 



for large momentum-space anisotropies and, in the meantime, focus on quantities that are 
independent of the freeze-out prescription. 
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(a) Wounded Nucleon 



(b) Binary Collision 



e x , 47ir|/S=1 
e p , 47ir|/S=1 
e x , 47it|/S=10 

e p , 47IT|/S=10 




e x , 47it|/S=1 
e p , 47IT|/S=1 
e x , 4jit|/S=10 
Ep, 4jit|/S=10 




FIG. 3. (Color online) Spatial and momentum eccentricities as a function of proper time for (a) a 
Glauber wounded-nucleon transverse profile and (b) a Glauber binary-collision transverse profile 
with b = 7 fm, A = T = 0.6 GeV, £ = 0, and u ±fi = at r = 0.25 fm/c. For the 4tt7]/S = 1 
run we used An = To = 0.6 GeV and for the 4nr]/S = 10 run we used An = To = 0.576 GeV 
for wounded-nucleon initial conditions and An = To = 0.584 for binary-collision initial conditions. 
These adjustments were made in order to guarantee the same final particle number. In all cases we 
used the centered-differences algorithm with a lattice size of 100 x 100, a lattice spacing of a = 0.4 
fm, and a RK4 temporal step size e = 0.01 fm/c. 



A. Smooth Initial Conditions 



We begin by presenting results using smooth initial conditions. For numerical tests of 
the various algorithms we refer the reader to App. [B} Therein we show scalings with lattice 
spacing, box size, and comparisons of the different algorithms employed for both smooth 
and fluctuating initial conditions. 

In order to demonstrate that AHYDRO reproduces known 2nd-order viscous hydrodynam- 
ics results, in Figs. [T] and [2] we compare the results of an aHydro run with results obtained 
using the latest version of the code of Romatschke and Luzum [12] . In Fig. [I] we assumed 
Airrj/S = 0.1 and in Fig. [2] we assumed Aixri/S = 10. In both cases we show the isotropic 
temperature profile, T iso = 7^ 1 / 4 (^)£i SO (A), in the left panel and the ratio of the y-component 
of the four velocity to the r-component in the right column. As can be seen from Fig. [T] there 
are only small differences at large radii in the case that the shear viscosity to entropy ratio 
is small. This demonstrates that our code reproduces 2nd-order viscous hydrodynamics in 
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T iso [GeV] at x = 0.50 fm/c T iso [GeV] at x = 1 .50 fm/c T iso [GeV] at x = 2.50 fm/c 
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FIG. 4. (Color online) Visualization of the isotropic temperature and pressure anisotropy at three 
different times after the nuclear impact. For these plots we assumed a non-central collision with 
6 = 7 fm, an isotropic Glauber wounded-nucleon profile, and a b = fm central temperature of 0.6 
GeV at 0.25 fm/c. For this plot we used a value of Attt]/S = 1 and a lattice size of 200 x 200 with 
a lattice spacing of a = 0.2 fm and a RK4 temporal step size of e = 0.01 fm/c. 

the limit of small rj/S. Fig. [2] shows the case of large shear viscosity to entropy ratio. In this 
case we see only small deviations in the temperature profiles and substantial differences in 
the flow profiles. We therefore expect the aHydro and 2nd-order viscous hydrodynamics 
frameworks to give different flow observables for large rj/S. We note that corrections near 
the edges are expected even for small values of rj/S and that the relative magnitude of the 
aHydro flow and the viscous hydrodynamics flow is to be expected: since aHydro gener- 
ates larger longitudinal pressure than viscous hydrodynamics one expects diminished radial 
flow. This pattern is also observed in simulations which use the lattice-boltzmann method 

In Fig. [3^i and [3]d we compare the spatial and transverse momentum-space eccentrici- 
ties as a function of proper time assuming two different values of the shear viscosity to en- 
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FIG. 5. (Color online) Visualization of the isotropic temperature and pressure anisotropy at three 
different times after the nuclear impact. For these plots we assumed a non-central collision with 
6 = 7 fm, an isotropic Glauber wounded-nucleon profile, and a b = fm central temperature of 0.6 
GeV at 0.25 fm/c. For this plot we used a value of Airrj/S = 10 and a lattice size of 200 x 200 with 
a lattice spacing of a = 0.2 fm and a RK4 temporal step size of e = 0.01 fm/c. 

tropy density ratio corresponding to typical strong-coupling {A-kt)/S = 1) and weak-coupling 
(47T7//iS = 10) values. In Fig. [3^ we used smooth Glauber wounded-nucleon initial con- 
ditions and in Fig. [3Jd we used smooth Glauber binary collision initial conditions. In both 
figures we assumed 6 = 7 fm, A = T = 0.6 GeV, £o = 0, and Uj_ : o = at r = 0.25 
fm/c and used the centered-differences algorithm with a lattice size of 100 x 100, a lattice 
spacing of a = 0.4 fm, and a temporal step size of e = 0.01 fm/c. In both cases RK4 with 
a temporal step size of e = 0.01 fm/c was used for the updates. As can be seen from these 
figures increasing the shear viscosity to entropy ratio by a factor of ten only decreases the 
momentum-space eccentricity e p at 5 fm/c by approximately 10% in both cases shown. We 
note, however, that the dynamical framework employed here, namely assuming that the lo- 
cal rest frame energy momentum tensor is azimuthally symmetric in momentum-space may 
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underestimate the full effect of the shear viscosity. 

In Figs. [4] and [5] we present visualizations in the form of colormaps with contours of the 
proper-time dependence of the isotropic temperature and the pressure anisotropy defined by 
the ratio of the longitudinal and transverse pressures. Fig. [4] shows the case of Airq/S = 1 
and Fig. [5] shows the case of Atiti/S = 10. In both cases we assumed a non-central collision 
with 6 = 7 fm, a Glauber wounded-nucleon profile, and a b = fm central temperature of 
A = T = 0.6 GeV at r = 0.25 fm/c. A lattice size of 200 x 200 with a lattice spacing 
of a = 0.2 fm and a RK4 temporal step size of e = 0.01 fm/c was used in both cases. As 
we can see from this figure the magnitude of the momentum-space anisotropics can be large 
in the center of the fireball and grows towards the edges. In Fig. [4] we see that assuming 
4irr]/S = 1 at r = 1.5 fm/c the center still has a 25% momentum-space anisotropy and 
assuming 4nr]/S = 10 (Fig. [5]) one finds approximately 85% momentum-space anisotropy at 
t = 1.5 fm/c. In fact, in the case of 4nr]/S =10 the system is highly anisotropic during 
the entire evolution. For such large shear viscosities the aHydro framework provides a 
dynamical framework which should be more reliable than the naive application of 2nd-order 
viscous hydrodynamics. 

In Fig. [6] we plot the momentum space eccentricity, e p , at the "freeze-out time" Tf as a 
function of the assumed impact parameter, b. For this figure we used a Glauber wounded- 
nucleon transverse profile with £0 = 0, and u±fi = at To = 0.25 fm/c assuming 47Ti]/S = 1 
and 4nri/S = 10 and a freeze-out temperature of Tf = 0.15 GeV. For the 4nr]/S = 1 run we 
used A = T = 0.6 GeV as the central temperature and for the 4irr]/S = 10 run we used 
A = T = 0.576 GeV in order to guarantee the same final particle number. We used the 
centered-differences algorithm with a lattice size of 200 x 200, a lattice spacing of a = 0.2 
fm, and a RK4 temporal step size e = 0.01 fm/c. The freeze-out time Tf was determined 
by finding the time at which the maximum isotropic temperature T lso dropped below the 
freeze-out temperature of Tf = 0.15 GeV. This figure shows that changing the assumed 
value of the shear viscosity to entropy ratio from one to ten only makes a difference of 8% 
in the peak value of the momentum-space ellipticity. We should note, as a caveat which we 
will emphasize again in the conclusions, that because we assume that the energy momentum 
tensor is azimuthally symmetric in the local rest frame this places us somewhere between 
a full blown viscous hydrodynamical calculation and ideal hydrodynamics. Therefore, firm 
conclusions will have to wait until results with a completely general ellipsoidal energy- 
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FIG. 6. (Color online) Momentum eccentricity at the freeze-out time as a function of impact for 
an isotropic Glauber wounded-nucleon transverse profile with £o = 0, and u± i o = at tq = 0.25 
fm/c assuming Tf = 0.15 GeV. For the Atttj/S = 1 run we used Ao = To = 0.6 GeV as the central 
temperature and for the inrj/S = 10 run we used Ao = To = 0.576 GeV in order to guarantee the 
same final particle number. We used the centered-differences algorithm with a lattice size of 200 
x 200, a lattice spacing of a = 0.2 fm, and a RK4 temporal step size e = 0.01 fm/c. 

momentum tensor are available. 

B. Fluctuating Initial Conditions 

For our fluctuating initial condition case we have implemented Monte-Carlo (MC) 
Glauber initial conditions [66]. At a given impact parameter b we statistically sample 
a Woods-Saxon distribution to determine the position of the nucleons in each colliding 
nuclei. We then compute the transverse distance between each pair of nucleons from nuclei 
A and B and assume that they collide if the transverse distance between the centers of the 
nucleons being compared is less than d = \Z~onnJtt. If a collision is deemed to have occurred 
a two dimensional gaussian with width cr = 0.46 fm is added to the energy density. We 
then adjust the overall scale to match the smooth Glauber model results. 

In Fig. [7] we present visualizations in the form of colormaps with contours of the proper- 
time dependence of the isotropic temperature and the pressure anisotropy defined by the 
ratio of the longitudinal and transverse pressures. In Fig. [7]we assumed a central collision b = 
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FIG. 7. (Color online) Visualization of the isotropic temperature and pressure anisotropy at three 
different times after the nuclear impact. For these plots we assumed a collision centrality of b = 7 
fm with a sampled Monte-Carlo Glauber wounded-nucleon profile and an isotropic temperature of 
T = 0.6 GeV at 0.25 fm/c. For this plot we used a value of Airrj/S = 1. We used a lattice size of 
200 x 200 with a lattice spacing of a = 0.2 fm and a RK4 temporal step size of e = 0.01 fm/c. 



7 fm with a sampled Monte-Carlo Glauber wounded-nucleon profile, an isotropic temperature 
of A = T = 0.6 GeV at r =0.25 fm/c, and Atttj/S = 1. We used a lattice size of 200 x 200 
with a lattice spacing of a = 0.2 fm and a RK4 temporal step size of e = 0.01 fm/c. As 
can be seen from this figure, fluctuations can induce large momentum-space anisotropies, 
particularly in regions where the initial temperature is lower and therefore the relaxation rate 
is smaller. In a 2nd-order viscous hydrodynamical approach one would have many "spots" 
with very large momentum-space anisotropies. Note that Fig. [7] shows the case Aixri/S = 1 
and we do not include a similar figure for the case of Airrj/S = 10; however, we note that 
similarly to the case of smooth initial conditions, for this large value of the shear viscosity 
to entropy ratio, one sees large persistent momentum-space anisotropies throughout the 
simulated region. 
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VI. CONCLUSIONS 



In this paper we studied the application of anisotropic hydrodynamics to the evolution 
of the matter created in relativistic heavy ion collisions. We began by specifying a tensor 
basis for the energy-momentum tensor which was applicable when the system is azimuthally 
symmetric such that one has energy density transverse pressure, and longitudinal pressure 
along the diagonal in the local rest frame. Microscopically we were able to demonstrate that 
if one assumes local momentum-space azimuthal symmetry, it suffices to introduce one scale 
A and an anisotropy parameter, £, which controls the transverse-longitudinal momentum- 
space anisotropy. 

We then used these results in the computation of moments of the Boltzmann equation. 
Using the zeroth and first moments of the Boltzmann equation we were able to determine 
dynamical equations for the plasma scale, A, anisotropy parameter, £, and the transverse 
flow components u x and u y . In order to solve the resulting partial differential equations 
we implemented three differencing schemes: centered differences, weighted LAX, and hy- 
brid Kurganov-Tadmor. The first method is suitable for smooth initial conditions whereas 
the second two are required when one considers event-by-event simulations. Based on our 
analysis and benchmarks we find the weighted LAX scheme to be faster than the hybrid 
Kurganov-Tadmor scheme with both giving the same results within controllable numerical 
errors. 

We showed through explicit solution of the resulting partial differential equations that the 
pressure components remain positive definite and that plasma momentum-space anisotropics 
grow larger as one approaches the transverse edge. In addition, we studied fluctuating initial 
conditions and demonstrated that fluctuations can result in regions of high momentum-space 
anisotropy in the center of the simulated matter. As a cross check we demonstrated that 
in the limit of small rj/S the solution of the AHYDRO dynamical equations reproduces 
results from publicly available 2nd-order viscous hydrodynamics codes. For smooth initial 
conditions we demonstrated that, subject to the assumption of momentum-space azimuthal 
symmetry in the local rest frame, one sees a relatively small variation of the final lab frame 
momentum-space eccentricity e p as rj/S is increased. Drawing quantitative conclusions from 
the results contained herein might be premature, however, since the impact of relaxing 
the assumption of azimuthal isotropy of the energy momentum tensor in the local rest 
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frame is unknown. Removing this assumption will result in what we will term "ellipsoidal" 
anisotropic hydrodynamics. Work in this direction is currently underway. 

We note in closing that there have been a number of authors studying the behavior of 
anisotropic plasmas in strongly coupled gauge theories [341 1361 167H71] . The aHydro frame- 
work agrees extremely well with existing 1st, 2nd, and 3rd order viscous hydrodynamical 
results which have been computed analytically for strongly-coupled M = 4 supersymmetric 
Yang-Mills |75J . It would be interesting to see if any of the results contained herein could be 
used in the context of strongly-coupled theories in order to develop useful phenomenological 
models. One open question first raised in Ref. [H] concerns whether or not the breaking of 
rotational symmetry in momentum-space requires the introduction of transverse and longitu- 
dinal transport coefficients. Mathematically this would seem to be the case in our formalism 
if one linearizes fluctuations around an anisotropic background. Such possibilities will be 
explored in the future. In the meantime, the progress made here opens up the possibility 
for phenomenological application to heavy ion observables such as collective flow, photon 
and dilepton production, quarkonium screening, jet energy loss, etc. in the presence of large 
momentum-space anisotropics. 
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Appendix A: Particle production in the (0+l)-dimensional case 

In this appendix we discuss the issue of particle production in 2nd-order viscous hydrody- 
namics vs anisotropic hydrodynamics. To begin we note that there are two limits in which 
one expects particle production to go to zero: (a) the limit of ideal hydrodynamics and (b) 
the free-streaming limit. For small but non- vanishing shear viscosity we expect there to be 
additional particles associated with dissipation; however, as the shear viscosity to entropy 
ratio increases we should see a maximum in the particle production since it will eventually 
have to go to zero in the free-streaming limit. In contrast, second-order viscous hydrody- 
namics predicts that the excess in particle production is a monotonically increasing function 
of the assumed value of rj/S. 




0.1 1 10 100 1000 10000 5 10 15 20 



47iT|/S 4jit|/S 

FIG. 8. (Color online) Total particle number at r = tj as a function of the assumed value of the 
shear viscosity to entropy ratio. For this figure we ignored transverse expansion making the system 
effectively (0+l)-dhnensional and we used initial values of Ao = 0.6 GeV and £o = at tq = 0.25 
fm/c. 

In order to demonstrate the difference quantitatively, in Fig. [8] we plot the quantity 
r/r n/riQ — 1 at r = Tf as a function of Anrj/S. We used a freeze out temperature of Tf = 150 
MeV to determine Tf. This quantity should be zero if there are no particles produced during 
the evolution. As can be seen from these plots our expectations are confirmed, namely 
that one sees a maximum in entropy production at large values of ^mr\jS with it returning 
to zero as Anri/S increases above this point. Concentrating on the zoomed plot in Fig. [8] 
one sees that for 47rr)/S = 10 2nd-order viscous hydrodynamics overestimates the entropy 
production by approximately 93%. We note that as the initial temperature is lowered, the 
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FIG. 9. (Color online) Spatial and momentum eccentricities as a function of proper time for a 
smooth Glauber wounded-nucleon transverse profile with 6 = 7 fm, Ao = To = 0.6 GeV, £o = 0, 
and uj_.o = at tq = 0.25 fm/c assuming 4ttt]/S = 1. In all three cases we used a RK4 temporal 
step size of e = 0.01 fm/c. 

excess particle production obtained from 2nd-order viscous hydrodynamics becomes larger. 
This will be important for phenomenology since one of the key constraints on rj/S stems 
from having to reduce the assumed initial temperature in order to compensate for dissipative 
particle/entropy production. 
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Appendix B: Numerical Tests 

In Fig. [9] we show the time evolution of the spatial and transverse momentum-space 
eccentricities as a function of proper time for a smooth Glauber wounded-nucleon transverse 
profile with 6 = 7 fm, A = T = 0.6 GeV, £0 = 0, and u± i o = at r = 0.25 fm/c assuming 
4m] /S = 1. In all three cases we used a RK4 temporal step size of e = 0.01 fm/c. In this 
figure we have used the central-differences algorithm without wLAX smoothing and compare 
the effect of varying the lattice spacing and lattice volume. As can be seen from this figure, 
the systematics are well under control in this case. Knowing that the centered-differences 
algorithm systematics are under control we can now compare with the hybrid Kurganov- 



Tadmor algorithm. In Fig. 10 we show such a comparison for the same conditions as shown 



in Fig. |9j As can be seen from this figure the naive centered-differences algorithm and the 
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FIG. 10. (Color online) Spatial and momentum eccentricities as a function of proper time for a 
smooth Glauber wounded-nucleon transverse profile with 6 = 7 fm, Ao = To = 0.6 GeV, £o = 0, 
and u±fl = at To = 0.25 fm/c assuming Airrj/S = 1. Here we compare the centered-differences 
and Hybrid Kurganov-Tadmor algorithms. 
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FIG. 11. (Color online) Spatial and momentum eccentricities as a function of proper time for a 
smooth Glauber wounded-nucleon transverse profile with 6 = 7 fm, Ao = To = 0.6 GeV, £o = 0, 
and u±fl = at to = 0.25 fm/c assuming Airrj/S = 1. Here we demonstrate the convergence of 
the wLAX algorithm with A = 0.05 to the result obtained without any spatial smoothing as one 
decreases the lattice spacing. In all cases RK4 with a temporal step size of e = 0.01 fm/c was used. 



e x , Centered Differences (300x300 a=0.2 fm) — i — 
e p , Centered Differences (300x300 a=0.2 fm) - - 
e x , Hybrid Kurganov-Tadmor (300x300 a=0.2 fm) - - -» - - - 




e x , Centered Differences (200x200 a=0.2 fm) — i — 
e p , Centered Differences (200x200 a=0.2 fm) - -k- - 
e x , Centered Differences + wLAX (200x200 a=0.2 fm) - - * - ■ - 



E p , Centered Differences + wLAX (200x200 a=0.2 fm) ....q.... 
e x , Centered Differences + wLAX (400x400 a=0.1 fm) 
E p , Centered Differences + wLAX (400x400 a=0.1 fm) - ■ o- - 
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FIG. 12. (Color online) Spatial and momentum eccentricities as a function of proper time for a 
sampled MC Glauber wounded-nucleon transverse profile with 6 = 7 fm, £o = 0, and u± t o = at 
ro = 0.25 fm/c assuming Attij/S = 1. Here we compare the hybrid Kurganov-Tadmor and wLAX 
algorithms. For the wLAX update we used RK4 with a temporal step size of e = 0.01 fm/c. 



hybrid Kurganov-Tadmor algorithm give results that are indistinguishable by eye. 



In Fig. 11 we present the spatial and momentum eccentricities as a function of proper 
time for a smooth Glauber wounded-nucleon transverse profile with 6 = 7 fm, Ao = To = 0.6 
GeV, £o = 0, and u± t o = at r = 0.25 fm/c assuming inrj/S = 1. In this plot we compare 
a run with the unsmeared centered-differences algorithm and the wLAX algorithm with two 
different lattice spacings. As can be seen from this figure the amount of numerical viscosity 
is small and can be reduced if one reduces the lattice spacing. 



To further illustrate the reliability of the wLAX algorithm in Fig. 12 we compare a single 
MC Glauber wounded-nucleon run using both the wLAX and Hybrid Kurganov-Tadmor 
algorithms. Both codes were initialized with the same sampled MC initial condition (a 
visualization of the evolution of this configuration is shown in Fig. [7]). As can be seen from 
this figure, wLAX and Hybrid Kurganov-Tadmor give virtually indistinguishable results. We 
point out in this context that the wLAX algorithm take much less time to complete a run 
giving it a significant advantage when one wants to sample many different configurations. 
Based on our benchmarks the wLAX algorithm is approximately ten times faster than the 
Hybrid Kurganov-Tadmor algorithm. 
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Appendix C: Boost invariant Id dynamics - The Bjorken solution 



In this section we briefly review what happens when the system is boost invariant, ho- 
mogeneous in the transverse directions, and has conserved particle number, i.e. Jo = 0. For 
this situation, it is convenient to switch to the comoving Milne coordinates defined as 

t = r cosh q , 

z = r sinh^ • (CI) 

In this coordinate system the metric g^ v = diag (1, —1, —1, — r 2 ). In addition, the local rest 
frame four-velocity simplifies to 

= (cosh c,0,0, sinh^) , (C2) 
such that u T — 1, u q — 0, and we have 

D = = d T , 

= duvP = - . (C3) 

T 



By applying the last two expressions to the zeroth moment of the Boltzmann equation (2.26) 
for an isotropic plasma we obtain 



which has a solution of the form 
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d T n = , (C4) 

r 



n(r) = n ^ . (C5) 
r 



If now we apply again the expressions given in Eq. (C3) to the first moment of the 



Boltzmann equation (Eq.2.32) one finds easily that 



E -+- V 

d T £ + — — = . (C6) 

T 

If the system has an ideal equation of state (EOS) then £ = 3V and one can further simplify 
this to 

dr£ = -\-, (C7) 



which has a solution 

r / 



<^idealgas — <^0 ( — ) • (C8) 
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If the system does not have an ideal EOS but instead has an equation of state corresponding 
to a constant speed of sound, i.e. dV/dS = c 2 , then it follows that V = c 2 s £ where we have 
fixed the constant by demanding that the pressure goes to zero when the energy density 
goes to zero. In this case one finds instead 

£ = £o(^) 1+< \ (G9) 

which reduces to the ideal case when c 2 s = 1/3. If the EOS has varying speed of sound then 
one can express V in terms of an integral of the speed of sound. Alternatively, one could 
calculate the pressure and energy density separately for e.g. an ideal massive Boltzmann 
gas [76J for which one finds 



£ = iVdof 2 

Z7T 

e /VT m 2 T 2 /; , 

V = N dof — — — K 2 



TJ 

2tt 2 ~"*Vt)' 

n=^, (CIO) 

and 

^"-^-H^T- (cu) 

Note that the thermodynamic relations above are consistent with Bjorken scaling for the 
number density, n/n = r /r, for all values of m in the case of isotropic hydrodynamics. 
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